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O Abstract 

\^ This is the second in a series of papers aimed at developing a practical time-domain method for self- 

iy-s force calculations in Kerr spacetime. The key elements of the method are (i) removal of a singular part of 

the perturbation field with a suitable analytic "puncture" based on the Detweiler- Whiting decomposition, 

'TT' (ii) decomposition of the perturbation equations in azimuthal (TO-)modes, taking advantage of the axial 

^i symmetry of the Kerr background, (iii) numerical evolution of the individual jn-modes in 2-|-l-dimensions 

I with a finite difference scheme, and (iv) reconstruction of the physical self-force from the mode sum. Here 

j^JT) we report an implementation of the method to compute the scalar-field self-force along circular equatorial 

'— ' geodesic orbits around a Kerr black hole. This constitutes a first time-domain computation of the self force 

^_ in Kerr geometry. Our time-domain code reproduces the results of a recent frequency-domain calculation by 

^ Warburton and Barack, but has the added advantage of being readily adaptable to include the back-reaction 

fSJ from the self force in a self-consistent manner. In a forthcoming paper — the third in the series — we apply 

^~H our method to the gravitational self- force (in the Lorenz gauge). 
O 
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I. INTRODUCTION AND OVERVIEW 

The context of this work is the ongoing program to apply the theory of self forces (SFs) in 
curved spacetime to the problem of a massive compact body orbiting a black hole [H [2]. This 
program is now reaching fruition. Recent highlights include a calculation of the gravitational SF 
along bound (but otherwise generic) geodesic orbits in Schwarzschild geometry j^, the quantitative 
determination of several gauge-invariant post-geodesic effects associated with the gravitational SF 
[IHS], and the first meaningful comparisons with post- Newtonian results 0^9]. Results from SF 
calculations are being used to inform the development of an Effective One Body (EOB) model of 
the inspiral [9l [10] , and even as a benchmark for fully- nonlinear numerical-relativistic simulations 
|llj . Despite this significant progress, gravitational SF study so far has been restricted to the 
relatively simple model where the background geometry is that of a non-rotating (Schwarzschild) 
black hole. In order for the program to become truly astrophysically relevant [12], the community 
must now come to address the problem of the gravitational SF in Kerr geometry. 

A first calculation of the SF in Kerr spacetime has in fact been presented recently |13t 114] . 
employing the scalar-field SF as a (relatively) simple toy model for the gravitational SF itself. This 
calculation — first performed for a circular orbit in the equatorial plane of the Kerr black hole |13j 
and later extended to eccentric orbits (still in the equatorial plane) |14j — is based on a frequency- 
domain treatment of the field equations: the (scalar field) perturbation equations are integrated 
numerically mode-by-mode in a Fourier-harmonic decomposition, and the physical scalar-field SF is 
then constructed using standard mode-sum regularization [15]. This approach is computationally 
efficient and relatively easy to implement (as it reduces the computational task to, essentially, 
the solution of ordinary differential equations). However, the method also has several drawbacks: 
it becomes rapidly less efficient with increasing eccentricity; it cannot be incorporated easily in a 
scheme that evolves the orbit self-consistently under the back- reaction effect of the SF; and it is not 
easily generalizable to the gravitational SF in Kerr. (An alternative frequency-domain approach, 
designed specifically to deal with the last of the above challenges, is being developed by Shah et 
al. [mdZ].) 

Time-domain methods offer a natural way around the above difficulties. Numerical codes based 
on time evolution are versatile and flexible, and can accommodate different orbital configurations 
with relatively little alteration. They can, in principle, be readily adapted to incorporate the 
back reaction from the SF on the orbit in a self-consistent manner. The metric perturbation 
equations in the Lorenz-gauge lend themselves most naturally to a time-domain treatment thanks 
to their hyperbolic nature. Lastly, time-domain codes formulated in 2+lD or 3+lD do not rely on 
separability of the perturbation equations and thus provide a natural route to the Kerr problem. 
For these (and other) reasons, the time-domain approach to SF calculations has been favoured by 
a number of workers in the field, and is under active development by several groups |18H24j . 

The obvious weakness of time-domain formulations lies in their computational costliness. The 
lack of separability of the perturbation equations into multipole harmonics on the Kerr background 
means that one has to deal with partial differential equations in 2+lD (or in full 3+lD). An 
added complexity is introduced by the presence of a pointlike source: the physical perturbation 
diverges along the worldline of a pointlike particle (the divergence is Coulomb-like in the 3+ ID 
case and logarithmic in the 2+lD case [25]), and a naive numerical implementation would fail. 
One first has to "regularize" the field equations in order to make them amenable to numerical 
treatment. A regularization scheme formulated in 2+lD was introduced in Refs. |25L I26j . and a 
closely related method is being developed by Vega and collaborators [HI |27l [28] within a 3+ ID 
framework. Both methods are based on the idea of a "puncture" function: Rather than evolving the 
physical (singular) field, one evolves a "punctured" version thereof, obtained by subtracting from 
the full field a certain (singular) puncture function, given analytically. In practice, this amounts 



to (numerically) solving the field equations with a certain "effective source" to obtain the residual, 
regularized field. If a suitable puncture function is chosen (see below), the physical SF can be 
constructed directly from the residual field. 

The 2+lD scheme of Refs. |25t [2B] takes advantage of the axial symmetry of the Kerr back- 
ground. The perturbation equations are separated into azimuthal (oc e*™'^) modes and solved 
mode-by-mode, using an ?7i-mode version of the puncture scheme. The SF is then reconstructed as 
a sum over m-mode contributions (see below for a more precise description) . In Ref . ^24j (hereafter 
'Paper I') we presented a first full implementation of this m-mode regularization scheme, as applied 
to the scalar-field SF on a Schwarzschild background. We used this simplified model as a conve- 
nient platform for development, and in order to illustrate the generic properties of the method in 
a relatively "clean" environment. 

In the current work we take an important step forward by applying our method to the Kerr 
case. We thus present the first time-domain SF calculation in Kerr geometry. As in paper I, we 
consider a scalar-field toy model, and specialize to circular orbits (in the equatorial plane of the 
Kerr black hole). This lays the groundwork for a gravitational SF implementation in Kerr, which 
we will present in the third paper of the series [29] . 

In the remainder of this introductory section we review the ?7i-mode regularization method in 
some detail, and give a short preview of the current work. 

A. ?Ti-mode regularization 

Below is a schematic description of the m-mode regularization method as applied to the scalar- 
field SF in an axisymmetric spacetime such as Kerr (the basic structure of the scheme remains 
unaltered in the gravitational case). We model the compact object as a point particle carrying 
a scalar charge q and moving along a geodesic of the background (e.g., Kerr) geometry. In our 
model, the scalar field $ sourced by the particle is taken to satisfy the minimally-coupled Klein 
Gordon equation, 

n$ = S, (1) 

where the covariant derivatives in D = V^V^ (as elsewhere in this work) are taken with respect 
to the background (e.g. Kerr) metric ga^, and the source S is that associated with the energy- 
momentum of the pointlike charge q [modelled by a delta-function distribution; cf. Eq. ( [I5| ) below]. 
In our notation, $ represents the "physical", retarded solution to Eq. (jTl) (which is, of course, 
divergent at the particle's location). Let us recall that the physical SF acting on the scalar charge 
can be derived by means of the Detweiler- Whiting formula ^30j 

F,%,, = qV''<^R = qV'^{'^-<^s). (2) 

Here ^s ("S field") is a specific singular solution of the inhomogeneous equation M, and the 
difference $/j = <I> — <I>5 ("R field") is a smooth, homogeneous solution of the sourceless field 
equation, O^r = 0. The formal construction of $5 is described in [30] . 

In most cases of interest, the singular field $5(x) cannot be evaluated in exact form. However, 
the form of ^^(x) near the particle's worldline can be studied analytically to inform a local approx- 
imation $77(2;) — the "puncture function" — which may be defined globally through some extension 
far from the worldline. The "order" of the puncture describes how well it approximates <I*5(x) 
near the particle. In Paper I we introduced the following convention: a puncture is said to be of 
"order n" , denoted $p , if the difference $p — ^s is a C"~^ function on the particle's worldline. 
In addition, we require $p = <I>s as well as V"<I>p = V^^s at the particle limit (noting that the 



first of these conditions makes sense only for n > 2, and the second only for n > 3). Since the R 
field $/j is smooth (C°°), the residual field associated with the n-th order puncture, 

ci>N^,j,_$W^^^_(^N_^^)^ (3) 

is a C"~^ function on the worldline. Hence, for example, a 2nd-order residual function <l>y is 

[31 

continuous but not differentiable, and a 3rd-order residual function $^ is both continuous and 

In] 

differentiable. For n > 3 we have V"$^ = V"<l>/j at the particle limit, which allows us to write 

F,:,f = gV"cI>t^^^l. (4) 

Thus, the problem of calculating the SF reduces to the problem of computing the residual 
field ^ij^ (for some n > 3) near the particle. This field satisfies a "regularized" version of the 
Klein-Gordon equation: 

n$gi = 5(x)-nci>^U52. (5) 

The effective source S^^{x) is an extended source that no longer involves a delta function; it is 
generally C"'~^ on the worldline. Hence, a 4th-order puncture gives rise to a continuous effective 
source, in general. Following the development of efficient computational tools for constructing 
high-order Green function expansions |28t [3H [32] , a 4th-order puncture for generic orbits in Kerr 
geometry is now at hand, and higher-order extensions are possible. 

The 3-l-lD treatment of Vega et al. [IHl ETJ [28] tackles the field equation ([s]) directly. In the 
2+ ID formulation, instead, we take advantage of the azimuthal symmetry of the background in 
order to reduce the dimensionality of the problem. Let (/9 be a suitable "azimuthal" coordinate 
associated with the axial symmetry of the background (below we will make a specific choice for (/? 
in the Kerr case; note that the Boyer-Lindquist coordinate </> does not suite our purpose as it is 
singular at the event horizon of the Kerr black hole). We formally decompose the residual field 
into azimuthal modes in the form 

oo 
m=— oo 

and similarly for <l>-p and S^s (we have omitted the labels [n] for brevity). This decomposition 
separates the field equation (15]) , with each of the ttz- modes satisfying an equation of the form 

n-$^ = 5,"^, (7) 



where D"^ is a certain m-dependent differential operator in 2-l-lD [cf. Eq. ( |25| ) below]. Given the 
modes <&^(j;) near the particle, the SF can be reconstructed as a sum over modal contributions, 

oo 

KeM= E ^s"el?, F^^ ^ qV^<^^e^^^) , (8) 

m=— oo 

evaluated at the particle. That the gradient indeed commutes with the modal sum (for n > 2) was 
shown in Ref. |26j. Ref. |26j also argued that the reconstruction formula ([8| applies, in fact, for 
any ra > 2, even though the original 3+ ID formula (pi) makes sense only for n > 3. Note that the 
individual modal contributions F^^ depend on the puncture order n; however, the total SF -Fg^jf 
is, of course, n- independent (for n > 2). 



In our method, the fields $^(x) are obtained by solving Eq. (|7|) numerically for each m using a 
finite-difference algorithm on a 2+ ID mesh. The boundary conditions for <1>^ depend, in principle, 
on the global extension of the puncture $-p. It is convenient to truncate the support of ^-p far 
from the particle, in which case the residual $7?, (and its modes $^) satisfy the usual "retarded" 
boundary conditions (outgoing radiation at infinity and purely ingoing radiation across the horizon; 
boundary conditions for any nonradiative modes are determined from regularity). This truncation 
is achieved in our scheme by introducing an auxiliary worldtube around the particle's worldline (in 



the 2-l-lD space), as we shall describe in Sec. IV A 



In Paper I we implemented the above "Tn-mode regularization scheme" to compute the SF 
acting on a scalar charge in a circular orbit around a Schwarzschild black hole (refraining from 
a multipole decomposition and working instead in 2-l-lD). We implemented our numerical code 
with punctures of second, third, and fourth orders. In all cases our numerical results were found 
to be in good agreement with those of previous calculations using frequency-domain methods or 
time-domain integration in 1-l-lD, thus reaffirming the validity of the reconstruction formula (JSJ) 
for n = 2-4. 

We used the simple Schwarzschild model in Paper I to explore some of the generic features 
of the m-mode method. The question of the convergence rate of the modal sum in Eq. (Is]) is a 
crucial one, since in actual numerical implementations one can only compute a finite number of 
modes (realistically, no more than ~ 20 modes using our current code). We made the following 
empirical observations (supported by heuristic arguments): (i) The dissipative piece of the SF has 
an exponentially convergent mode sum; (ii) the conservative piece of the SF has a mode-sum that 
converges with an inverse-power law: the individual modal contributions in the mode sum fall off 
at large m as m~^ for even n and as m"""*"^ for odd n (where n is the order of the puncture used). 
Thus, the large- ttt, modal contribution to the conservative piece of F^^^ is oc m"^ with a 3rd-order 
puncture, and oc m~^ with either a 4th or a 5th-order puncture. We concluded that a 4th-order 
puncture implementation offers an optimal balance between simplicity (of the analytic puncture 
formulation) and efficiency (of the computational algorithm) . We expect this conclusion to carry 
over to the Kerr case. 

In our implementation (as in any other time-domain implementation) the initial data for the 
time evolution are, of course, unknown. Instead, we begin the evolution with arbitrarily selected 
initial data, and rely on the property of homogeneous perturbations in black hole geometries to 
completely "dissipate away" at late time ( "no hair theorem" ) . We let the evolution proceed until 
the solution relaxes to a stationary state at a sufficient level. It is important to understand the 
rate of relaxation, since this determines the necessary evolution time and thus has important effect 
on the computational burden. As discussed in Paper I, the theoretical expectation is for the modal 
contributions to relax with an TTi-dependent inverse-power law in time t. Modes with higher m 
relax faster, and the slowest relaxation is exhibited by the m = mode: oc t~^ for the field modes 
$^ and oc t~^ for the force modes F^^. This behavior was confirmed numerically in Paper I; 
we expect it, too, to carry over to the Kerr case. The issue points to an important advantage of 
the 2-l-lD formulation (as compared with the full 3-l-lD treatment): Since different ?TT,-modes are 
evolved separately, we have the flexibility of being able to adjust the evolution time as a function 
of m to achieve computational saving. 

B. This paper 

Here we extend the analysis of Paper I to the Kerr case, focusing again on circular motion (in 
the equatorial plane). Our 2-|-lD treatment makes such a generalization straightforward in prin- 
ciple. However, there are several nontrivial technical details that require careful consideration and 



extra development. First, the analytic expressions involved in the puncture formulation become 
considerably more complicated and require commensurably greater care in their numerical imple- 
mentation. Second, the azimuthal coordinate used to define the m-mode decomposition must be 
chosen judiciously; the standard Boyer-Lindquist coordinate suits for purpose in Schwarzschild ge- 
ometry but not in the Kerr case. The third point is most crucial: We find that the finite-difference 
scheme employed in Paper I to evolve the field equation luh becomes unstable in the Kerr case and 
must be replaced with an alternative method. We have experimented with several alternatives, 
and will report here our findings. We observe that some of these methods, used previously in the 
literature to study vacuum perturbations on Kerr, develop instabilities at large values of m, which 
makes them unsuitable for us. We develop and implement a finite-difference scheme (a variant 
of the "method of lines") that appears to perform well even for large m values. Our scheme is 
second-order convergent (in the grid spacing), and we implement it here with a puncture of 4th 
order (n = 4). 

The remainder of this paper is structured as follows. In Sec.lTIlwe describe our physical setup and 



review the formalism of scalar-field perturbations on a Kerr background. In Sec. Ill we formulate 
the m-mode regularization scheme as applied to the problem at hand, and prescribe 2nd, 3rd 
and 4th-order puncture functions with corresponding effective sources. Section |IV| describes our 
numerical method, including the worldtube technique and the finite-difference scheme. We also 
report on some of our less-successful experiments with other schemes. In Section |V] we present 
some numerical results along with convergence tests, and explore the late-time relaxation of the 
numerical solutions and the convergence properties of the ?TT--mode sum. We demonstrate that our 
code successfully reproduces the results of Warburton-Barack's frequency domain calculation |13j . 



In the concluding section. Sec. |VI[ we discuss future directions, including the application of our 
method to the gravitational SF. 

Throughout this work we use metric signature — 1-++ and set G = c= 1. 

II. PRELIMINARIES 
A. Setup and notation 

In Boyer-Lindquist coordinates (t, r, 9, (p), the line element ds"^ = g^ydx^dx'^ of a Kerr spacetime 
with mass M and spin parameter a is given by 

,2 A 2Mr\ ,2 4aMrsin2 6l, , , /O^ , 2 1 ,ni ( 1 2 IMra^ &\y? e\ . 2 ^,,2 
ds^ = -\\ 2-) '^^ 2 dtd(\>^j-dr^^p^dQ^^\r^ ^c? ^ ^ ) sm^ 6I#^ 

where 

p^ = r^^c?co^e and A = r^-2Mr + a^. (9) 

The spacetime admits two Killing vectors, ^^r. = S^ and ^f^^ = ^a,, as well as a Killing tensor Q^"^. 
The event horizon is the hypersurface r = r^, where r_|_ is the larger of the two roots of A(r), 
given by r-t = M it \/M'^ — a^. The ergosphere (inside which static observers do not exist) extends 
between the event horizon and the static limit surface Vgi = M ^ \J M"^ — a? cos^ Q. Stationary 
observers just outside the event horizon have angular velocity {d(j)/dt) given by 

which may be interpreted as the angular velocity of the horizon. Timelike geodesies of the Kerr 
geometry can be parametrized by the three constants of motion 

£ = -?('t)'"M) ^z = i1^)U^, Q = Q^'^Uf.Uy, (11) 



referred to as the (specific) energy, azimuthal angular momentum and Carter constant, respectively. 
In the last expressions u^ = g^yu'^ ■, where u'^ is the tangent four- velocity along the geodesic, and, 
as elsewhere in this work, we use the background (Kerr) metric to raise and lower indices. 

In this work we consider the scalar-field SF on a pointlike scalar charge moving on a circular 
geodesic orbit in the equatorial plane of the central Kerr black hole. The motion may be either 
prograde (i.e., in the same sense as the rotation of the black hole) or retrograde. We shall adopt 
the convention that u'^ (hence also Cz) is always positive, and distinguish between prograde and 
retrograde orbits by the sign of a: a > for the former; a < for the latter. 

For geodesies in the equatorial plane {9 = 7r/2) one finds Q = 0. If the orbit is both equatorial 
and circular, the other two constants of motion are given by 

l-lv"^ + av^ ^ 1 - 2av^ + a?v^ ,^ „. 

c = —;= , Cz = rnv , =, (12) 

where rg is the Boyer-Lindquist radius of the orbit, v = ^ MJtq and a = a/M. The particle's 
angular frequency is given by 

^ = ^ = !f! = I ns) 

dt M* M{l + ~av^y ^ ' 

B. The scalar-field equation in 2-|-lD 

The scalar field is assumed to be governed by the sourced (minimally coupled) Klein-Gordon 
equation, 

n$(x) = ^gd^ [./^gg^^'dMx)] = S{x), (14) 

where g = —p'^siv?6 is the Kerr metric determinant, and g^'^ is the contravariant form of the 
metric. The source S{x) is taken to have a delta-function support on the geodesic worldline: 

/oo 
[-g[x)]-'^/^5^[x-x{T')]dT'. (15) 

-oo 

Here q is the scalar charge of the particle, r is its proper time, and x = x{t) [shorthand for 
x" = x"(r)] parameterizes the particle's geodesic orbit. In our setup we have f = ro(=const) and 
9 = 7r/2, as well as ^ = fit (taking ^ = at t = without loss of generality). 



We now wish to decompose the field equation ( 14 ) into azimuthal modes. We must be mindful 
of the fact that the Boyer-Lindquist azimuthal coordinate (f) is pathological at the event horizon. 
We cure this by making the standard coordinate change (j) ^ (p defined through 

dip = d(l) + -7-dr, (16) 

which, choosing the constant of integration so that ip coincides with (p as r ^ oo, gives 



¥?((/), r) = (f)-\ In 



r_ 






(17) 



Since both <l> and S are periodic in ip, they admit m-mode decompositions (formally, Fourier 
expansions) of the form 

oo oo 

$= ^ ^mgim^^ S= Y^ S-^e"^^. (18) 

m=— oo m=— oo 



The ?TT,-modes ^"^ may be obtained via the inversion formula 

1 



#™(t,r,^) 



27r 



e-'"''^^{t,r,e,ip)d^, 



and similarly for S"^. For the latter we obtain, performing the integration, 



S'' 



rlu 



-S{r-ro)5{9-TT/2)e- 



-irrnp(t) 



(19) 



(20) 



where (p{t) = (p((j)(t),rQ). Here we used the notation (j){t) = (j){T{t)), where the function r(t) is 
obtained by formally inverting the (monotonically increasing) function i(r). 

The above decomposition separates Eq. (14) into a set of 2+lD equations for the individual 
(complex- valued) ?7z-modes $"^(t, r, 9). To write these equations in a form more amenable to numer- 
ical treatment via the Method of Lines (see Sec. IV B below), we introduce two minor modifications: 
we transform to the new variable 



and we replace dr — )• dr, , where the Kerr tortoise coordinate r* is defined (as usual) by 



or, specifying the constant of integration 

2M 

r^: = r -\ 



dr 



r+ In 



r — rj^ 



2M 



r_ In 



r — r- 



2M 



r^ — r 
This brings the separated wave equation into the final, working form 

□?\ 

where the m-mode scalar wave operator reads 

^™ 9^ AiamMr d 

■^ " 9t2 + S2 dt 
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■S2 



5^ + cote 



de 



□m^m ^ ^m^ 








- reads 






(r2 + a2)2 d^ 


"2iamr(r2 + a2) -2a'^A' 
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the m-mode source is given by 



and we have introduced 



am 



rt\p^ 



5" 



(21) 



(22) 



(23) 



(24) 



(25) 



(26) 



r^ + a^)^-a^Asin^6l. 



(27) 



III. PUNCTURE FORMULATION 



The retarded solutions of the m-mode field equation (24) diverge on the particle's worldline for 



all m. [The divergence is logarithmic in the spatial distance to the particle, and its leading-order 



form does not depend on m — cf. Eq. (23) of Ref. [25].] In our approach we do not tackle Eq. (24) 
directly, but rather consider a "punctured" version thereof, given by 

r-\m,j,[n]m _ (,ln]m ^„Q^ 

Here the variables ^^ ™ are defined as the m,- modes of the nth-order residual field 

xi/N=r(c^_ci,W)=^,j,W^ (29) 



i\n\m 



41fr = -^(^-°4"')- (30) 



with S'ljjgg being the m-modes of the associated nth-order effective source 

^[n] __r 

With a suitable puncture ^^~ , the physical SF is then constructed via 

oo 

i^s"eif= E ^s"el?, Ff,I? = gV" (r-iM';^['^le*'"^) , (31) 

m=—oo 

evaluated at x(t). To complete the formulation one need only prescribe a suitable puncture $^ , 
which can then be used to construct (analytically) the effective source S^j,^ via Eq. (30). 

In] 

In what follows we construct puncture functions <l>p of 2nd, 3rd and 4th orders, based a local 
expansion of the Detweiler-Whiting S field. Similar punctures have been developed in previous 
work [261 EH |33] , but we give them here (and in Appendix |A]) in explicit ready-to- use coordinate 
form, specialized to the problem at hand of circular equatorial orbits in Kerr. We then also describe 
the procedure for constructing the m modes <l>p and S^^^, and give explicit analytic expressions 
for n = 2. 



A. Punctures through fourth order 

Our first task is to obtain a local expansion for the Detweiler-Whiting S field <&s'(x) in the 
vicinity of a point on the worldline, x. Let us presume that x and x are connected by a unique 
spacelike geodesic, so that we may introduce the Synge bi-scalar a{x,x) [HIS!, which is equal to 
one half of the squared geodesic distance along the geodesic connecting x and x. Furthermore, we 
presume that there exists a sufficiently-regular coordinate chart which covers a region of spacetime 
encompassing x, x and the spacelike geodesic. In this chart, points x and x take coordinates x" 
and x", and the vector tangent to the worldline at x is denoted by u°'. Now we let ct^ denote the 
covariant derivative of a taken with respect to x, and we introduce the orthogonal and parallel 
projections 

^2 = (^a/3 ^ -a-^^-^-^^ ^ _ -a-^^ (33) 

where g'^^ is the background metric evaluated at x. Note that e has the interpretation of the spatial 
geodesic distance from the field point x to the worldline, i.e., the length of the spatial geodesic 
segment connecting x to the worldline and normal to it |47] . 



At this point we introduce the coordinate difference 5x, with coordinates 



(5x" 



X 



(33) 



Our aim below is to obtain a coordinate expansion for the Detweiler- Whiting S field, ^six"), in 
powers of dx°'. When Sx"' is sufficiently small, ^s can be shown to admit the local expansion 
[281 [33] 



^six) 



1 + 



a[x, x) 



+ 



/3(x,x) 



+ 0{6x^) 



where 



a = - (?7^ - e^) Rai3^s{x)u"v?a'^a\ 



(34) 



(35) 



and 



/3 



24 



7?(7?^ - 3e^)V^i?„/3^5(x)n"n^n^CT^a' - (r/^ - e')V ^Ro,^^5{x)u'^u^^'^^^ 



<^r,'y?rP?r^Frl^ 



(36) 



Here Rap-yS and 'Vfj^Rap-yS are the background (Kerr) Riemann tensor and its covariant derivative 
(both evaluated here at x), and ct" = g'^^aj^. Note the bi-scalars a and (3 are smooth (C°°) 
functions of x and x (and hence 5x), with a = 0{5x'^) and /3 = 0{5x^) near the worldline. Hence, 
in Eq. (34) the terms a/e^ and j3/e^ are, respectively, C^ and C^ functions of 5x at the worldline. 
Therefore, defining the "truncated" S fields 



$ 



[2] 



9/e, 



$ 



[3] 



/(£-' + 



ae 



$ 



[4] 



q{e + ae 



Pe- 



rn 



we find that the differences $5 — ^^ are C"~^ at the worldline. 

The local behavior of the functions ^^ ' ' near the worldline is consistent with that of punctures 
$p' ' (correspondingly). However, in this form ^g (x) cannot be used in a practical scheme for 
two reasons: (i) (Ja [and therefore $^'(x)] is not known in explicit analytic form and (ii) aa [and 

therefore $^ (x)] is only defined if the field point x is in the vicinity of the worldline point x. To 
construct globally-defined analytically-given punctures we take two steps: (i) We expand e^, a and 
/3 in powers of 6x, to obtain a local coordinate expansion of ^g up to a suitable order, (ii) We 
promote the local expansion to a global function by choosing an appropriate analytic extension. 
To take step (i) , we use the method of Ref . [3l] to obtain a Taylor expansion for aa in powers 



of dx (with fixed x). We use the Taylor expansion in Eqs. (32), (35) and (36) to obtain expansions 



for the bi-scalars e^, a and /?. Let us use the notation S(„), «(„) and /3(„) to denote the Taylor 
expansion of the corresponding quantities up to and including the 0{5x") term. 

To take step (ii), we allow the point on the worldline x to become a function of the field point 
X, i.e., x(x), in such a way that x and x(x) are connected by a unique spatial geodesic as x 
approaches the worldline. A simple way to achieve this is to take t = t, where from here on we 
specialize to Boyer-Lindquist coordinates. Namely, we relate a field point x° = (t, r, 9, (j)) to a 
worldline point x° = (t, tq, '7r/2, ilt). Then the coordinate differences are simply 5t = Q, 5r = r — vq 
and 59 = 9 — 11/2. We note that (j) — Vtt\s not periodic, and it may be large even when the points are 
close together. To rectify this we instead choose 54> = 2sin~"'^ {sin[((/) — r2t)/2]} where the principal 
values — 7r/2 < sin~^ x < 7r/2 are assumed. 

After taking steps (i) and (ii), we obtain truncated Taylor series coordinate expansions S(„), 
«(„) and /3(-„) in place of the geometric quantities e^, a and /3. Explicitly, e^ = S(„) + 0{5x'^~^^)^ 
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where 



5(„) 



m2 Y1 s,,km\69yi6^)\ (38) 



2<i+j + k<n 



and likewise for a and /3. Here, for convenience, we have introduced the dimensionless radial 
distance 6r = {r — ro)/M, and Sijk are dimensionless coefficients depending only on tq and a. In 
Appendix IA] we give explicit expressions (for the specific case of circular equatorial geodesies on 
Kerr spacetime) for all coefficients Sijk necessary for computing S(„) through n = 5. We also give 
explicit expressions for «(„) and /3(-„) through n = 5. 

With the quantities 8(3^4^5), a{4,5) and /3(5) all at hand in explicit analytic form, 2nd, 3rd and 
4th-order punctures can be prescribed as follows: 

*?' = -^. (39) 

5(3) 
5(4) ^(2) 



*? = Wn^ + ^ + ^l- m 




Note that in these expressions we have relaxed the assumption that 6x°' is small in any sense, and 
allow X to be any field point in the spacetime outside the black hole. (In practice we shall only need 



X to lie within a certain finite- width worldtube surrounding the worldline — see Sec. IV A below.) 
Since we work in a particular coordinate system (i.e., Boyer-Lindquist), this procedure amounts to 
choosing a particular global extension of $5 off the worldline. The above construction guarantees 

\n\ 

that our globally prescribed punctures <I>p approximate the S field $5 at the required level near 
the particle's worldline, i.e., that $p = $s up to a C^~'^ function. 

B. TO-mode decomposition of the puncture field 

We now describe the construction of the puncture m-modes $p for n = 2 and n = 4. We 
skip the n = 3 case, which will not be considered further in this work. 

1. Second- order puncture 

The method for obtaining the ?7z-mode decomposition of a 2nd-order puncture field of the form 



[21 

(39) was described in Paper I. The first step is to recast $p as a periodic function of 6<j), without 



any loss of order n. For n = 2, this may be achieved by making the replacement 

<5(/.2 ^ 2(1 - cos 6(1)) = 6^^ + {5(1)^) . (42) 

Note that the puncture function involves only even powers of 6cl), and so it suffices to replace 6(j)'^ . 



Next, following the steps of Paper I [and using Eq. (19)] we write 

1 r ■ [91 g-imA</) 

27r./_ ^ ^ 2tt 



^[2]m ^ ^ / g^im^^g]^^ ^ ^ / e-™^$gl#, (43) 
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where we have introduced 



A(/) = y?- 



In 



r+ 






(44) 



[recall Eq. ([17|]. The integrand in Eq. (|43| has the form (g/M)e~^"'<^ [^ + 25(1 - cos 5(/>)]"^/^ 
with 



A = 5200(5^ + 3020^6* + SsooSr-^ + 5i2o6f5e , B = S002 + Sl02<5r, 



(45) 



where the coefficients Sjjfc are those given in Appendix [Al This integral can be expressed in 
terms of the complete elliptic integrals of the first and second kinds, respectively K{k) = 



Jq (l — k"^ sin^ x) dx and E{k) = /J^ (l — k"^ s\v? x) dx, giving 

-Jm(™+A0) 



l-p 



where 



27rMSV2 
[A/(45)]i/2^ 
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(46) 
(47) 



The quantities p^{p) and p'^{p) are polynomials (of order m) in p^; these polynomials are tabulated 
for m = 0-5 in Table I of Ref. [25] and can be readily obtained for higher values of m. 



2. Fourth- order puncture 



.w 



To recast $p as a periodic function, without diminishing the order n = 4 of the expansion, we 
made the replacement 

5 8 _ 1 



cos 0(/) + - cos 25(j) = 6(j)'' + 0(6(j)^ 

2 3 6 



(48) 



Note that the 4th-order puncture, like the 2nd-order puncture, is symmetric under 6(j) — )• —S(f), and 
hence it features only Scp"^ and S^i"^ terms. In our implementation we chose to replace 6(1)^ with 



the square of Eq. (48). An alternative choice would be to replace 6(j)^ with a linear function of 



harmonics cos{k6(j)) (where A; = 0, 1, 2) via 



6 — 8 cos 5(j) + 2 cos 25(j) 



+ 0((5/ 



(49) 



The m-mode decomposition proceeds as in Eq. (43). In this case, however, the form of the 



fourth-order puncture given in Eq. (41) is such that the integration cannot be readily written in 



terms of elliptic integrals. Instead, we proceed by direct numerical integration of Eq. (43). 



Effective source 



The effective source 5'^gg is found by applying the wave operator D to the puncture [see Eq. ( 30 )] , 
recalling that the dependence of <I>p on t and is only through the combination </> — J7t (for circular 
equatorial orbits). In practice, this implies that 

d ^ d d d ^^^^ 



dt 



-n 



d{6(t)) ' d(l) 5((50) ■ 



i[n]m 



The m-mode quantities S\j,^g are then obtained via an inversion formula analogous to Eq. (19). 



In the 2nd-order case, the integrals are readily computed analytically in terms of elliptic Integra 
(see Appendix pi). In the 4th-order case, we found it to be more straightforward to compute the 
integrals via direct numerical integration. 
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IV. NUMERICAL METHOD 

In this section we describe the key features of our numerical implementation. These features are 
identical for any puncture order, so we shall omit here the order identifier '[n]' from our expressions, 
for brevity. 

A. Worldtube scheme and the numerical domain 

Our puncture function $-p and effective source Sx^eS are defined through a naive extension of 
a local expansion, and their behavior far from the particle may well prove problematic. Indeed, 



the 4th-order puncture in the form given in Eq. (41 ) generally diverges at r — )■ oo, and so does the 



effective source associated with this puncture. To handle this behaviour, we construct (following 
Ref. [25] and Paper I) an auxiliary "worldtube" 7" in the 2+ ID domain to enclose the worldline. 
Outside the worldtube, we evolve the finite-difference (FD) version of the homogeneous wave 



equation (24) for the full field modes ^'". Inside the worldtube, we evolve the FD version of the 



sourced wave equation (28) for the residual field modes ^^. Across the boundary of the auxiliary 
tube, dT, we convert between ^™ and ^^ using the analytically-given value of the puncture field. 
To summarise, 

n™*™ = 0, outside T, 

n™^/™ = 5^,g, inside r, (51) 

The full field modes ^™, of course, satisfy the usual physical boundary conditions at the asymptotic 
domains (see below). 

Our FD scheme is formulated on a fixed uniform grid based on t,r*,0 coordinates; see Fig. [T} 
The grid spacings in the corresponding coordinate directions are denoted At, /\r^, A6. In the 
2-t-lD domain, the particle's trajectory traces a straight line at 9 = 7r/2 and r^ = r=Ko[= r*(ro)]. 
For convenience, we lay the grid so that the particle's trajectory crosses through a row of grid 
points (of fixed r^:,6). 

On this grid, we use a worldtube T of fixed coordinate widths {F^, , Tg} centered at the worldline. 
Consider an arbitrary grid point with coordinates (t, r^, 6). If |r^,— r=Ko| < ^r^f^ and \9—7r/2\ < Tg/2 
then the point is said to lie within the worldtube; otherwise it lies outside. For convenience, we 
choose Tr^ and Tg to be integer multiples of the grid spacings Ar* and A6, respectively. 

Within the FD scheme described below, the value of the field at a given grid point (say, point 
"o") is computed based on the values at several neighboring grid points, obtained in previous 
steps. If all these grid points lie either inside or outside T, then the FD scheme is implemented 
straightforwardly. However, if point o lies close enough to dT that some of these points are "in" and 
others are "out" , we make the following adjustment. If the point o is "out" , then we first demote all 
"in" points to "out" points using ^"^ = ^^-|-r$p, before applying the FD formula. If, conversely, 
the point o is "in", then we promote all "out" points using ^^ = ^"^ — r^p, before applying 
the FD procedure. In this way, we make use of the known value of ^^ in order to communicate 
between the two numerical variables ^"^ and ^^ across the boundary of the worldtube. 

B. Finite-difference scheme: Method of Lines 

We describe our FD scheme as applied either inside or outside T- For concreteness, we refer 



specifically to the full field equation (24), but note that the same procedure is applied inside T, 
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^^ and 5^ 



qm 



with the replacements ^"^ 

For our purpose, we rewrite the m-mode field equation (1241) as two coupled first-order equations 
in the (complex-valued) variables ^"^ and 11'" = d^'^/dt. Explicitly, 



^l 



nr 



n*^ 



AiamMr 



n'" + 



(r2 + a2)2 



S^ 



ij/'" + 



2iamr{r'^ + a^ 



2a^A 



A ^^ A cot _^ _ 



A 

S2 



m 



sm 



+ 



2M 



-— ^ 






2iam 



^' 



s:^ 



* ) 



where ^™, ^J?, ^^ etc. denote partial derivatives. This form is suitable for the application of the 
method of lines, which we adopt here. (The method is reviewed, e.g., in Ref. [35] and in Sec. 4.2 
ofRef. |M]-) 

Let ^njfc = '^"^{tn,r^:j,6k) denote the field value at the grid point corresponding to the coordi- 
nate values 



tn = nAt, 



r*o + jAr*, and 9k = kA9. 



(52) 



To write Eqs. (52) in an FD form, we replace ^"^ — )• ^njk, and replace all spatial derivatives (in 



r* and 6 coordinates) with finite-difference approximations of second-order accuracy: 












^ 



n{j+l)k 



^ 



n{j-l)k 



2(Ar,) 

^n(i+l)fc - 2*„jfc + ^n{j-l)k 



{Any 



^ 



nJik+1) 



^ 



nj(k-l) 



2(Ae) 

^nj(fc+l) - 2^njA: + ^nj(k-l) 

(Key ■ 



(53) 
(54) 

(55) 
(56) 



To evolve the equations in time (i.e., to step forward to obtain {^ (^n+i)jkj^{n+i)jk} from {^njk^ 
Hnj/fc}) we applied a fourth-order Runge-Kutta step. To illustrate the method in general, consider 
a coupled set of first-order equations of the form 



dtu = F{t,u), 



(57) 



where u and F are vectors. We denote the set of unknowns u at time t„ as u„ (in our case, 
u„ = {^njk, n^jfc} for all values of j, k). To obtain an approximation Un+i at time tn+i = in + At 



we apply 

u„+i = u„ + At (6iki + b2k2 + bsks + biki) , 

which combines information from a sequence of intermediate estimates, 

ki = F(t„,,u„), 

k2 = F(t„ + C2At,u„ + Ata2iki), 

ks = F(t„ + C3At,u„ + At(a3iki +a32k2)), 

k4 = F {tn + c^At, u„ + At (a4iki + a42k2 + a43k3)) . 



(58) 



(59) 
(60) 
(61) 
(62) 



We used a standard version of the 4th-order Runge-Kutta method, in which the coefficients are 



0.21 



a32 



2, 043 — 1; 0-31 



041 



042 



0, 6i 



h C2 



C3 



and 
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FIG. 1: Grid for finite difference metfiod. Tlie left plot shows the 2+lD grid with linear spacings (Ai, Ar*, 
Ad). The evolution proceeds from initial data on the surface t = (shaded). The right plot illustrates a 
grid cell for the method-of-lines finite difference method (with the 9 direction suppressed for clarity). To 
obtain the value at point (j, n + 1), we first use the values at the 5 grid points {j — 2, n), . . . , (j + 2, n) to 
obtain estimates on an intermediate time slice (unfilled circles); cf. Eq. (58). The domain of dependence of 
the grid point at (j, n + 1) is shown with (densely) dotted lines. 



C4 = 1. Note the function F(-) couples between {^njk,'n^njk} and {^n{j±i){k±i),'^n{j±i){k±i)}; see 
the right panel of Fig. [T| which illustrates the method. 

The above scheme is second-order convergent in the spatial resolution and fourth-order con- 
vergent in the time resolution. We opted for a 4th-order convergence in time scheme because the 
method of lines with a 2nd-order Runge-Kutta time step is known to be unconditionally unstable 
[36], and (through experimentation) we found that a 3rd-order scheme requires an impractically 



small time step to remain stable. We comment more on stability in Sec. IV D below. 



C. Initial and boundary conditions 



We specify a "zero" initial condition on the initial surface at t = 0, i.e., ^l 



Ojk 



n-, for all 



j, k. Of course, this is not a solution of the sourced field equations; however (as discussed in Paper 
I and Ref. [37j) initial "junk radiation" is found to radiate away at late times. We monitor the level 
of this transient radiation on the worldline (by measuring the deviation from stationarity), and 
discard the early part of the evolution, where the magnitude of this radiation is high. The residual 
"contamination" (due to imperfect initial conditions) diminishes with time, but may remain a 
significant source of numerical error. We discuss this further in Sec. |VA 
The physical boundary conditions to be applied at the poles (0 = and 6 



^ 



m=0/ 



G{0,7r})=0, 



To implement these conditions, we simply set ^^ 
we extrapolate to obtain 



vr) are [25], 

^™^°(0G{O,7r}) = O. (63) 

: n'" = at the poles for m 7^ 0, and for m = 



^^"=0(0 = 0) 



- r4^'"=0( 



vr 



[4^ 



m=0/ 



A0) - ^^"=0(0 = 2A0)] + O {{/\ef) , (64) 

AO) - ^'"=0(0 = vr - 2A0)] + O ((A0)^) , (65) 

because ^"^ is 



TT 



and likewise for H™. Note that the error term is O((A0)^), rather than 0{{A9)^ 
an even function oi 9 at 6 = 0. 

Conditions are also required at the boundaries of the grid in r^,, denoted nmin and nmax- A 
correct but technically demanding approach is to impose "absorbing" boundary conditions, i.e.. 
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purely ingoing radiation at the horizon and purely outgoing radiation at spatial infinity. We found 
it difficult, however, to control sufficiently well the level of initial spurious radiation bouncing off 
the spatial boundaries and contaminating the late-time solution. Instead, we use a simple "zero- 
derivative" condition, ^™ = 11™ = at the boundaries r^min and r^nnax, and merely ensure that 
the grid is of sufficient span in r^, that the initial burst of spurious radiation does not have time to 
reflect from the boundary and re-encounter the worldline. This approach is less computationally 
economical but has the advantage of simplicity. 

D. Numerical stability of method 

Our FD strategy here differs from that of Paper I, where we used a semi-characteristic grid in 
coordinates u = t — r^,, v = t + r^, and 6 (in Schwarzschild) . Let us refer to this as the ^^uvO^^ 
method. It is less natural to apply the uv9 method in the Kerr case, where the presence of a first 
time derivative in the field equation complicates the scheme considerably. In fact, we were not able 
to obtain sufficiently stable evolutions based on the uvO grid and thus abandoned the approach in 
favour of the method of lines. 

In the Schwarzschild case we found that the uvO method was subject to a numerical instability 
(manifested in an exponential growth of high-frequency noise) if the angular spacing /\9 was set 
to be too small. More specifically, the method was found to be unstable if |M] 



A0 _ 1 
~h 



< -max ('r-^A^/^) ^/l + m?/A, (66) 



where h = Aw = /\v is the characteristic grid spacing, m is the azimuthal mode number, and in 
the Schwarzschild case A = r^ — 2Mr. The instability was seen to arise first near the poles, and 
we showed that it could be mitigated by moving the polar boundaries inward (see Sec. III.C.5 of 



Paper I). The version of the method of lines described in Sec. IV B suffers from a similar angular 
instability. Empirically, we find that an instability arises if 

^ < max (r-'^^^'A y^l + m2/4, (67) 

for At = Ar*. Again, moving the polar boundaries inward (particularly for the large-m, modes) 
mitigates this problem, as described in Paper I. 

Explicit FD methods may also suffer from "radial" instabilities if the ratio v = At/Ar^, is taken 
too large. We chose to implement the method of lines with a 4th-order Runge-Kutta time step 
primarily for its known good stability properties ^36j. We were able to use the ratio v = 1 without 
problems. Below we comment on some experiments with other methods. 

E. Alternative finite-difTerence methods 

We have experimented with several alternative FD methods for evolving the scalar field on the 
Kerr spacetime in 2-l-lD, some of which have already been used in the literature. We report here 
briefiy on our findings. 

Krivan et al. \6'6\ I39| (see also Sundararajan et al. [^) rewrote Teukolsky's master scalar-field 
equations in 2-|-lD as a coupled set of first-order equations in variables ^™ and H"*, where 

fjium 2 _|_ 2 a^T;m 



16 



They applied a two-step Lax-Wendroff method |41) to obtain a second-order-accurate algorithm to 
evolv e forw ard in time. In our experiments with this method, we found that an angular instability 
(Sec. |IVD[ ) arises if A9/Ar^ < TTma^{r-'^A^/'^)y^l + m'^/8 (for a = and Ai = Ar*), and 
becomes rapidly more severe with increasing \a\. Note that at large m this necessitates fixing the 



ratio A0/Ar* at larger values than in our method of lines [compare with Eq. (67)]. To achieve a 
given resolution in the angular direction therefore requires smaller values of Ar^, (and At), which 
increases the computational cost. 

In addition, we tried a version of Krivan et a/.'s approach using a leapfrog method [41j in place 
of the two-step Lax-Wendroff. With this approach, we found that, at least in the Schwarzschild 



limit, the angular stability condition (67) applies, just as in our method of lines. However, we also 



found that in the Kerr case (a 7^ 0), relatively small values of u were required to obtain stable runs. 

We have also experimented with replacing the 4th-order Runge-Kutta step in the method of 
lines with either (i) a 3rd-order Runge-Kutta step, or (ii) an iterative Crank-Nicholson step |41] . 
In either case, the stability properties were found to be more restrictive than with the 4th-order 
Runge-Kutta method. 

Our numerical experiments have been restricted to coordinate grids in t,r^:,9. Other authors 
have highlighted benefits of working with alternative coordinate foliations of black hole spacetimes 
[I2III3]. For example, Zenginoglu 01j has discussed how, by making a suitable replacement of the 
time coordinate, t = t + a{r), one may use the method of lines on constant-i hypersurfaces that 
are asymptotically null and intersect the event horizon and future null infinity. Since the fields ^"^ 
tend towards constants on such surfaces, one may use a compactification step. Racz and Toth 
recently presented their own implementation of a similar idea for the Kerr spacetime. 



Simulations, data extraction and numerical error 



For a given mode m we evolve the initial data (Sec. IV C) according to the finite difference 
scheme (Sec. IV B) with boundary conditions (63)-(65) on a 2+ ID grid (Fig. [I| in the range 
< t < imax- The numerical solution depends on "physical" parameters, rg, a and m, and a set of 
"numerical" parameters, {num.} = {At, Ar=K, A9,Tr,,T0,tmax, ■ ■ •}• We refer to a simulation for a 
given m, tq, a with a unique set of {num.} as a "run". 

In Sec. HID of Paper I, we described the procedure for computing the SF from a set of runs; 
let us recap the steps here. First, we ensure that tmax is sufficiently large that the field inside the 
worldtube has settled into a stationary state (see Sec. M below). Next, we read off the following 
estimates for the modal contributions to the r and (p components of the SF (for m > 0) : 



Frit) 
mt) 



Ao 
-1 



5.*^^ 



r^^i!"^ 



'n 



{t,ro,7r/2,nt), 



2mro"^Im[^^(t,ro,7r/2)e 



imQt~\ 



(69) 
(70) 



where Aq ^ A(r = ro), and where we have defined the (real) "total ?Ti--mode contribution" via 



trmp 



[y^(t,r, 0), m = 0. 

Here all quantities are evaluated at grid points on the worldline at late times t < tmaxi and the 
derivative with respect to r^, is found via central differencing. In principle, the SF components are 
then found via the mode sums 



pself 



m=0 



m 
r ; 



F: 



self 



E^^' 



(72) 



m=0 
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along with F^'^^^ = —QFf^^^ and, by symmetry, Fg*^^^ = 0. In practice, the modal values (69)-(71) 
extracted from a particular run depend somewhat upon the set of numerical parameters |iium.|. 
In other words, the modal values contain numerical error, which we define as the difference between 
the numerical solution and the (unknown) exact solution. 

Several key sources of numerical error were discussed in Sec. HIE of Paper I. To recap, we have: 
(i) discretization error, associated with the use of a finite grid spacing At, Ar^,, A9; (ii) worldtube 
error, which varies depending on the dimensions of the worldtube, r^^,re; (iii) relaxation-time 
error, associated with nonstationary residual spurious radiation; and (iv) source cancellation error, 
arising from small inaccuracies in the numerical evaluation of the effective source near the worldline 
due to delicate mutual cancellation of terms. There arise additional errors in computing the mode 
sums: (v) ?7i-mode summation error, due to imposing a finite large-m, cut-off m-max, and fitting 
the remainder with an approximate model; and (vi) mode cancellation error, which can arise from 
delicate mutual cancellation of individual m-mode contributions leading to a high relative error in 
the mode sum. 

In Paper I we described in some detail how to monitor, control and (to an extent) mitigate 
these sources of error. Our treatment is very similar in the Kerr case. In Sec. |V] we will illustrate 
the effect of the various sources of error in our particular Kerr implementation. 

G. Computational resource 

Numerical error can be reduced simply by applying greater computing resources to the problem. 
For example, errors (i) and (ii) diminish as we increase the resolution (i.e., decrease At, Ar*, A6), 
and error (iii) diminishes if we run the simulation for longer (i.e., increase tmax)- However, the 
runtime scales as t'^g^-^/{AtAr^:A9) and so in practice we soon encounter a practical limit. 

Luckily, our numerical algorithm lends itself easily to parallelization, with each run (for given 
m, a ro, {num.}) being assigned to an independent thread. To facilitate this parallel computation 
we make use of the Iridis 3 HPC resource. The SF is computed by post-processing the results of 
multiple runs, as described in the next section. Typically, to compute F^'^^^ and Ff^^^ for a given 
ro and a, we computed the modes ?7i = 0, . . . , 19 up to tmax = 300M at four different resolutions, 
ArjM = 1/8, 1/16, 1/24, 1/32 [with fixed ratios At/Ar* = 1 and AO/Ar^ = 7r/(6M)]. Each of 
the 80 runs is assigned to a separate node. Some additional resource is devoted to the m = mode, 
to ensure a fuller late-time relaxation of this slowest-relaxing mode (see below); typically, we run 
this mode up to tmax = 500M. 

The longest runtime for one complete simulation in this work (i.e., one computation of the SF 
with given tq, a) was ~ 24 hours. 

V. RESULTS AND ANALYSIS 

In this section we present a selection of results from our numerical simulations, and we discuss 
the challenge of minimizing numerical errors. We compare the results with those obtained in |13] . 
and find good agreement to within the estimated numerical error. 

A. m-modes 

Let us start by considering a typical run, i.e., a single simulation with given ro, a, and m, and 
with a particular set of numerical parameters. We can visualize the run by extracting data along 
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FIG. 2: Field modes (to = 0, 1, 2, 5) on a constant time slice (at t — t,nax) for circular orbits at tq = lOM, for 
a range of Kerr parameters (ajM — —0.9, —0.5, 0, 0.5, 0.9). The left plots show field modes at fixed d — 7r/2 
and the right plots show field modes at fixed r = rg. Inside the worldtube we show both the residual field 
^^ (forming the 'trough') and the full retarded field 'I'™, which is divergent on the worldline. We note that 
the rotation rate ajM has only a subtle effect on the field profile. 



three slices: (i) t = imax) & = 7r/2, i.e., in the equatorial plane, (ii) t = tmax,?"* = ^*o> i-e-, from 
pole to pole, crossing the worldline, and (iii) r+ = r*o, = vr/2, i.e., along the worldline. 

Fig. ^ shows typical ?7i-inode contributions to the field along the constant-t slices (i) and (ii), 
for an orbit at ro = lOM and a range of Kerr parameters a. The worldtube is visible as a central 
'trough'; inside the tube, we show both the residual ^^ and the full field ^™ (which diverges 
logarithmically as r — )• tq, — )• vr/2). These plots are similar to those for the Schwarzschild 
implementation (see Fig. 4 in Paper I). We note that the effect of black hole rotation upon the 
field mode profiles is quite subtle, although it has a more profound effect on the SF. 



Figure 



shows plots of ^^, F^ and FV^ as functions of t on the worldline [i.e., on slice (iii) 



for runs with tq = lOM, a = 0.5M and modes m = 0, 2, 4 and 6. After an initial burst of 



junk radiation (due to imperfect initial conditions. Sec. IV C), the modal quantities settle towards 
steady-state values. Visible in the figures are two types of transients: Initially, there are regular 
high-frequency oscillations (for m ^ 0) which may be identified as quasi-normal ringing (indeed, the 
ringing frequency is proportional to m as expected, and the exponential decay rate of the ringing 
seems roughly independent of m, also as expected). At later times the modes exhibit a second type 
of transient behavior: a power-law decay with an m-dependent exponent. In Paper I (Sec. IVA5 
with, e.g.. Fig. 10) we explored this power-low behavior in some detail, and demonstrated that, by 
fitting the decay of the field with an asymptotic model, we can extrapolate to i — >■ oo to extract a 
steady-state value. We implement the same method here. 
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FIG. 3: Typical time series data for the modes of the field ^-ji (left column), the radial SF Fr (center 
column) and the angular SF F^ (right column) for a range of m modes (top to bottom: m = 0,2,4,6). 
These data are taken from simulations with radius rg = lOM and Kerr parameter a = 0.5M, with numerical 
parameters At = Ar* = M/32, A9 = 7r/192 and F^. = 2.5M, Tg = 7r/4. Note that F™=° is identically 



B. Convergence with grid resolution 

For simplicity, let us fix the ratios /\r^/ /\t and A6/At and vary only At. We will refer to 
some quantity X{At) as being "fcth-order convergent" if X{At) converges to some finite value as 
At — )• in the following way: 



X(At) = X(At = 0) + 0{At''). 



(73) 



To test the rate of convergence k of some quantity X{At) [e.g., the residual field modes ^^(At) 
or the SF modes F^>{At)], we may take the ratios of the differences between the results of runs 
with different resolutions. A standard test is to examine the quantity 



x(At) 



X(2At) - X{At) 



X{At) -X{At/2)' 
If the method is A;th-order convergent then xi'^t) will approach 2 in the limit At — t- 0. 



(74) 
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4.17 


3.89 




6 


3.96 


3.94 




9 


3.82 


3.99 


0.5 


3 


4.17 


3.94 
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3.97 


3.85 




9 


3.82 


3.74 


0.9 


3 


4.19 


4.00 




6 


4.00 


3.83 




9 


3.77 


4.14 



TABLE I: Sample convergence tests for rg = 6Af , a — 0.5M, showing numerical values for x as defined in 



Eq. (74) for the median grid spacings Ar* = At = Af/16 and A^ = 7r/96, with t,„ax = 300M. 



In our implementation of the method of lines we use 2nd-order-accurate finite differencing on 
the spatial slices. Hence, we expect to find that ^"^ is 2nd-order convergent. In addition, we expect 
the SF quantities F^ and FV^ also to be 2nd-order convergent. In Table |I| we give a few sample 
values of x from our simulations, for a median resolution of At = M/16, and ratios l\r^/ /\t = 1 
and A6/At = tt/(6M). We find that these values fall in the range 3.7 < x < 4.3 at this resolution, 
suggesting a second-order convergence of our FD scheme, as expected. We verified that, as the 
resolution is increased (i.e.. At, Ar*, A^ are decreased) the results of the convergence test improve 
(i.e., X becomes closer to 4). We find that, in general, the convergence deteriorates somewhat at 
large values of ?n, > 15. 

Once we have confirmed that the code is 2nd-order convergent, we may apply a Richardson 
extrapolation ("Richardson's deferred approach to the limit" [41j) to reduce the discretization 
error. Figure |4] illustrates the idea. The left plot shows numerical results for the field on the 
worldline, as a function of simulation time t, for various resolutions 1/16 < x = At/M < 1/32 
(and a = 0). The right plot shows the values at t = 300M as a function of resolution, and compares 
with a similar data set obtained using the uv9 method (which is only valid in the a = case). By 
fitting the data set with a simple polynomial model, i.e. 

^-(X) = ^^(x = 0) + C2X^ + C3X^ (75) 

we may extrapolate to infinite resolution {x = 0). FigureH^shows that the extrapolated values from 
the two data sets are in excellent agreement. The difference between the extrapolated values from 
two data sets (for example, runs with different worldtube dimensions) provides a rough estimate 
of the error in the fitting procedure. 

C. Mode sums 

To compute the total SF, we first compute the values for a range of m modes (typically m = 
0, . . . , 19). As discussed in the Sec. Ill for circular orbits the modes of the angular component 
F^ (which is purely dissipative) are expected to converge exponentially fast, F7^ ~ exp(— Am). 
The upper left plot of Fig. [51 showing —FV^ as a function of m on a semilog scale, confirms this 
expectation. 

By contrast, the modes of the radial component Fr (which is purely conservative) are expected 
to exhibit power- law convergence, according to the rule F^ ~ m~^, where C = 2 for 2nd- and 3rd- 
order punctures, and C = 4 for 4th-order punctures. In Paper I, we obtained results for punctures 
of 2nd, 3rd and 4th order, and demonstrated the expected convergence rates (see Figs. 13 and 14 
in Paper I). In the Kerr case, we have only implemented the 4th-order puncture. The upper right 
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FIG. 4: Convergence of the FD scheme and Richardson extrapolation. The left plot shows the m — 2 mode 
of the residual field on the worldline, Vl/^, for a = and a circular orbit at tq = 6A/, as a function of time. 
The lines show the results from simulations with various grid resolutions, x — 1/16, 1/20, 1/24, 1/28 and 
1/32, where x = At/A/ = Ar*/A/ = GAd/n. The right plot shows the values at tmax = 300M as a function 
of X. The red circles on the upper line show the values taken from the simulations in the left plot. The 
blue squares on the lower line show the values from the uv6 method (which only works for a = 0). The 
two data sets have been fitted with cubic polynomial models, ^^ — ca + C2x'^ + c^x^ (red and blue lines). 
Reassuringly, the values Cq extrapolated from the two data sets (which correspond to the zero grid spacing 
limit) are found to be in very good agreement. 



plot of Fig. ^ shows the magnitude of the modes, \F^\, as a function of m, on a log-log scale, for 
ro = lOM and a range of a. This graphically confirms that F™ ~ ?ti-~^ in the large-?7z limit, as 
expected. The bottom plot in Fig. p] shows the rescaled variable, m^F^, which tends to a constant 
in the large-m limit. We see that the value of the Kerr parameter a does not affect the convergence 
rate. 



D. Self-force results 



Let us now present sample numerical results from our 4th-order puncture implementation. In 
Tables II and III we give the radial and angular components of the SF, F^'^^^ and F^^^, for circular 
orbits at three radii, tq = 6M, lOM and nsco- Here nsco denotes the radius of the innermost stable 
circular orbit, which is a root of the equation 



6Mr 



0. 



(76) 



We give SF results for a range of Kerr parameters, a/M = —0.9, —0.7, —0.5, 0, 0.5, 0.7 and 0.9. 

The lower value in each entry of Tables O and III was obtained by Warburton and Barack |13| 
using Z-mode regularization in the frequency domain. The upper figure in each row was obtained 
via the ?7z-mode scheme described in the previous sections, with figures in brackets indicating the 
estimated error bar on the last quoted decimal. (The results from [13J are said to be accurate to 
all digits quoted.) Unsurprisingly (for our circular-orbit implementation), the time-domain figures 
are significantly less accurate than the frequency-domain results. 

Each (upper) value in Tables [ll| and III was obtained by post-processing the results from multiple 
runs. First, we used runs of various resolutions to extrapolate to zero grid spacing (see Sec. VB). 
For the modes m = (with the largest relaxation error) we further fitted the late-time data with 



a power-law relaxation model. To account for the large-?7i tail (Sec. VC), i.e., the modes m > 19, 
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FIG. 5: Convergence of modal contributions as a function of m, for a 4th-ordcr puncture implementation 
(ro = lOAf, tniax = 300M). The upper left plot shows (minus) the m modes of the angular (dissipative) 
component i^^, on a semilog scale. We observe exponential convergence with m. The upper right plot shows 
the m modes of the radial (conservative) component Fr, on a log- log scale. The dotted reference line is 
oc m~^. The lower plot shows the same data but rescaled by the expected power-law for a fourth-order 



puncture, i.e., m'^F^. We show a range of values of the Kerr parameter, a/M 
we observe that rotation does not affect the convergence rate. 



-0.9,-0.5,0,0.5,0.9, and 



we fitted the modes m = 12, . . . , 19 with a three-term model Am + Bm~^ -\- Cm~^. Estimates of 
the residual errors that remain after performing these steps were obtained. To estimate the total 
numerical error we combined the error estimates from discretization, relaxation, and tail-fitting 
steps in quadrature. 

We find that the results of the m-mode method are in good agreement with the frequency domain 
results of Ref. [13 j , and the difference between the data sets lies broadly within the estimated 
error bars of our results. The error budget is dominated by discretization and relaxation errors, 
with tail- fitting error subdominant in our 4th-order puncture implementation. We hope that, 
in future, implementation of a truly 4th-order FD scheme (using 4th-order spatial differencing) 
will diminish the discretization error to the level where late-time relaxation error is left as the 
dominant error source. Methods for dealing with the latter (such as hyperboloidal slicing [33], or 
multigrid refinement [Mj) are also ready for deployment. Hence it may be possible to arrive at a 
situation where tail- fitting is once again the dominant error source, in turn motivating the use of 
yet-higher-order punctures. 

It is worth noting that, though F^ is positive for all circular orbits on Schwarzschild (and all 
retrograde orbits on Kerr), it is negative for some prograde orbits on Kerr. This phenomenon was 
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4.102(1) 5 
4.100712 ^" 


1.1077(2) 
1.107625 


xlO"'' 
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-0.5Af 
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a — 





1.6771(2) 
1.677283 


xlO-'' 


1.379(1) 5 
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X lO-'"^ 
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-9.528095 


xlO-5 


-1.0913(9) 5 
-1.091819 ^^^ 


-1.0886(4) 
-1.088457 


xlO-3 
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+0.9M 


-1.6458(5) 
-1.645525 


xlO-4 


-1.767(1) 5 
-1.768232 ""^^ 


-1.1344(9) 
-1.133673 


xlO-2 



TABLE II: Numerical results for the radial component F^'^^\ for a range of Kerr rotation parameters a and 
orbital radii tq. Here nsco is the orbital radius of the innermost stable circular orbit (ISCO), which is a 
function of a and takes the values nsco/M = 8.717352,8.142965,7.554585,6.0,4.233003,3.393128,2.320883 
for a/M — —0.9, —0.7, . . . , 0.9, respectively. The table compares the results of our implementation of m- 
mode regularization in the time domain using a 4th-order puncture (upper entries) , with the results from the 
Z-mode regularization method in the frequency domain from Warburton and Barack }13lll4j (lower entries). 
Parenthetical figures in the upper values indicate the estimated error bar on the last quoted decimals. All 
digits shown in the lower values are significant. Entries left empty correspond to orbits below the ISCO. 
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TABLE III: Same as in Table [nj for the angular component F^"^^. 

clearly illustrated in Fig. 5 of [T3] , and is evident in the results of Table M The condition F^'^^^ = 
describes a curve in the ro-a plane. We note that the relative accuracy of the m-mode sum result 
diminishes for ro and a near this curve, as the individual m mode contributions conspire to cancel 
(though the sum may be zero, the individual m- modes are not). This means that small relative 
errors in individual m modes become magnified into large relative errors in the mode sum. In 
Paper I, we described this numerical problem, so-called "mode cancellation error" , in the context 
of orbits at large radii. 
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VI. OUTLOOK 

This paper reports on a first time-domain calculation of the self force in Kerr spacetime. Our 
purpose here was to develop a computational infrastructure (puncture functions, numerical methods 
and a working code) suitable for the Kerr problem, and illustrate the efficacy of our m-mode 
regularization approach. 

There are a number of obvious ways in which the efficiency of our numerical method may be 
improved. One may implement a higher-order-convergent FD scheme, and/or use a higher order 
(n > 4) puncture. The former improvement would lead to more rapid convergence with respect to 
the FD grid spacing, while the latter would accelerate the convergence of the m-mode sum. We 
recall our result from Paper I, according to which the individual modal contributions to the SF 
fall off at large m as m~^ for even order n, but only as 7n~^^^ for odd order n. Hence, it would 
be desirable to go directly to n = 6. The derivation of such a high-order puncture is possible by 



extending the formulation of Sec. Ill A to higher order using the methods developed in Ref. j31j : 
such an extension would be straightforward in principle, even if tedious in practice. Moving to 
a higher-order-convergent FD scheme would necessitate either (i) dealing more carefully with the 
nonsmoothness of the residual field along the worldline, or (ii) using a higher-order puncture in 
order to improve the differentiability of this field. Other possible ways to improve the efficiency 
of the numerical method include (i) implementing a mesh-refinement algorithm such as the one 
described in Paper I (or the more systematic algorithm developed by Thornburg |201 [2T] ) , or (ii) 
improving the treatment of boundary conditions, e.g., via the hyperbolic slicing method of Ref. 



Several natural extensions of this work are possible. First, one may consider more generic or- 
bits: eccentric, inclined, or even unbound; our time-domain framework has the advantage that it 
requires only a little adaptation to accommodate different classes of orbits. The puncture formula- 



tion described in Sec. Ill A is applicable for any type of motion in Kerr spacetime, with the specific 



orbital details entering only through the explicit values of the coefficients of S(„), a(„) and P(ny The 



latter are easily determined for any given orbit based on the covariant expressions (32), (35) and 



(36). The generalization of our FD algorithm to cope with generic orbits is also rather straightfor- 
ward: With a 4th-order puncture implementation, the residual field is twice differentiable at the 
particle (cf. Table I of Paper I), and — better still — its individual ttz- modes, which are the numerical 
variables, are thrice differentiable there (and smooth anywhere else). This means that grid points 
in the vicinity of the worldline would not require any special treatment within the FD scheme 
even when the orbits are non-circular (as long as one is content with a 2nd-order-convergent FD 
scheme, as we do in the current work). This stands in sharp contrast to the situation in the 1-l-lD 
framework of Ref. [3], where the lack of differentiability of the numerical variables made the FD 
scheme significantly more complicated in the case of eccentric orbits, and the generalization from 
circular orbits required much further development. 

Once a generic-orbit scalar-field SF code is at hand, one may attempt another important ex- 
tension: the inclusion of the back-reaction effect of the SF on the orbital motion, to compute the 
evolving orbit in a self-consistent manner. The time-domain framework allows, in principle, to 
correct the motion "in real time" as the evolution proceeds. However, it is a non-trivial task to 
formulate an algorithm for incorporating the SF information in a manner that is both physically 
robust and computationally efficient. The scalar-field SF problem offers a convenient environment 
for test and development, as it is devoid of gauge-related complexities. A self-consistent evolution 
has not been attempted so far, neither for the gravitational SF nor for the scalar-field SF. 

Of course, the most important extension of our code is to the gravitational case. We have 
already begun work to apply the strategy of ?TT,-mode regularization to the gravitational SF. The 
generalization of the puncture formulation to the gravitational case is very straightforward, and 
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gravitational-field punctures akin to $p are already available through 4th-order. Our numerical 
strategy, following [3l [18], is to work directly with the components of the Lorenz-gauge metric 
perturbation as numerical variables. The linearized Einstein equations are decomposed into m- 
modes, yielding (for each m) a set of 10 coupled partial differential equations in 2+lD for the 
various coordinate components of the m-decomposed metric. The key element of our method 
involves a judicious use of the (m-decomposed) Lorenz gauge conditions in order to reformulate 
the set of 2+lD perturbation equations so that they are amenable to numerical time-evolution and 
so that gauge-condition violations are automatically suppressed during the evolution. Care should 
also be taken to ensure that suitable boundary conditions are used for the various components 
at the poles. However, the numerical algorithm is otherwise similar to the one developed in the 
current work. We have already been able to successfully implement this approach in the test case 
of circular orbits on Schwarzschild, reproducing some of the gravitational SF results of Ref. |18j . 

However, using our method, we have thus far only been able to solve for the modes m >2. The 
two lowest modes, m = and m = 1, which contain nonradiative monopole and dipole pieces, do 
not evolve stably using our numerical scheme: for these two modes we find our numerical solutions 
to be growing linearly in time (even though gauge violations remain small). This situation is 
similar to the one described in the 1+lD framework of Refs. [31 [18], where no stable evolution 
scheme has been found for the monopole (/ = 0) and dipole (/ = 1) modes. Instead, Refs. |3l [18] 
had to rely on frequency-domain methods to obtain the monopole and dipole contributions. Our 
ambition is to find a method for calculating the two problematic modes (m = 0, 1, in our case) 
using time-evolution in 2+lD, relieving us from the need to resort to frequency-domain methods, 
the latter being useful only for bound orbits with sufficiently small eccentricity. We have made 
significant progress towards this end, and will report our findings in a forthcoming third paper of 
this series, to be focused on the gravitational SF problem in Kerr. 
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Appendix A: Coefficients for puncture functions through fourth order 



Here we give explicit expressions for the coefficients Sijk appearing in Eq. (38), specialized to 
circular equatorial geodesic orbits of Boyer-Lindquist radius tq about a Kerr black hole of mass 
M and spin parameter a. We give all coefficients necessary for computing S(„) through n = 5. 
We also give explicit expressions for the functions a(„) and /3(„) through n = 5. This constitutes 



sufficient input for constructing the 2nd, 3rd and 4th-order punctures through Eqs. (39)-(41). For 
compactness, we give the coefficients in terms of the dimensionless quantities 

fo = ro/M, a = a/M, (Al) 

in terms of which we have 

v = r-^^\ Ao = M-^A{ro) = rl-2fo + a^ (A2) 

as well as 

u'^^Mn'^ = ^^ ' u'= ;(^o + «^) . (A3) 

^ov ^0 — 3 + 2av yfo — 3 + 2av 
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Let us begin with the coefficients Sijk of the e^-expansion (38). The non-zero coefficients at 
0{Sx'^) are given by 



S002 = Ao('u 



,~,t^2 



S020 = ^'o 



S200 



at 0{6x^) by 



S102 = (^0 - l)(^i ) , S120 = ro, 



S300 






fo(d2 - ro) 



(A4) 



An2 



(A5) 



at 0{6x^) by 



S220 



fo — 3a 
6A0 



7^2 



S040 



12 L 



3a -fo(fo-2) 



S400 



2a^fo(l - 6fo) + 3a4 + fo^(8fo - 1) 



12An3 



1\q{u ) , \[ / \ —1 2 2/ N 

soo4 = ,^^ ^ [av - roj a(ro - 5)f + 2a + ro (ro + 1, 



S022 



S202 



6fo 
6t;2Ao 



12fo^ 
a^(9fo + 11) + 253t;-i(fo + l)(3fo - 8) + d'^roi^fo^ + lO^o^ - 3fo + 2) 

+2aT;-5(3fo' - H^o + 2) + ro\ro - l)(3ro - 2)] , 
ro - 3 + 25t;) [6Ao(a + v"^)^ + 5^(9fo + 7) + 2a^v~^(3fo^ - 7fo - 4) 



+d'fo{fo' - l)(3fo + 5) - Udv-'{ro - 1) - fo^(ro - 1) 



(A6) 



and, finahy, at 0{5x^), they are given by 



S500 



I,, -, a2(5-6fo)+ro^ 
S140 - ^(1 - ro), S320 ^^Ao^ ' 

5^(4 - 9fo) + a2fo(12fo^ - 5fo + 3) - fo^(6ro^ - 2fo + 1) 



12An4 
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{u'i'f 



12fo' 



Sl04 



a^(9fo + 29) - 2a^t;-^(5fo + 16) + d^roiSro - 17fo^ + 12) 
+2afo^t;-^(ll - 6fo) - 2fo^(3fo^ + 2fo - 3)] , 

fo\fo^ + 4fo - 9) + at;-9(fo' - 8fo + 9) + d''ro^i8fo^ + 3fo - 21) 
a^v"^(fo + 3)(6fo - 11) + S^ro(3fo^ + 15fo - 10) - a^v"^(3fo + 19) + 6a^ 






12fo 



f„4 



S302 = — ^^ (ro - 3 + 2~av) fo'={fo - 1) - 4Aofo + a^fo[4Ao(3fo - 4fo + 5) - 6fo 

4 , in~ 3 , c~ 2 ex; 1 ;;6/o~ , ^\ , r,;c5„,-l/~ , ^\ ~4 



-WOT x9fo^ + 6fo^ - 6fo] - a^(3fo + 7) + 2a57;-Hro + 4) - a^(23fo^ - ro^ - 26f 
+4a^7;~^(2Ao - 3fo^ + lOfo^ - 8fo) + 2at;"^[8Ao + fo(13fo^ - 22fo + 6 - 12Ao) 

For the 3rd-order puncture [Eq. (40)] we additionahy require the function au\, and for the 



(A7) 
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4th-order puncture [Eq. (41)] we require a(5) and /3(5)- These are compactly given in the form 

Aofo^Sr^ + Ao^h^Se^ + Ao^fo[a\ro + 2) + fo^]S^^} 



M^iu't'? 



(4) 






+ Ao^ro{d + v 



-3\2r ,2 



(A8) 



a 



(5) 



6a7^ 



-{Aofo^5f2 ^ Ao2fo^<5^2 + Ao^ro[d\fo + 2) + fo=^]5</.2 + [fo'(a' - ro)6f^ 

+ AoVo^-J^^ + Ao^(fo^ - a^)(5(A^]5f||Aofo^(2at>-^ - 35^ + 3fo - 2fo^)(5r2 

AoVo^(9d2 - lOat;"^ 



3\2r ,2 



+ AoVo^(-4a«-i + 3a2 + fo^)^^^ + Ao^fo(a + w" 

+ 4fo^ - 5fo)5e^ + Ao2(a2fo(3ro^ - 2fo + 5) + 2ai;-5(fo - 1) - 3a^ + ro\4fo - 7) 

+ ^0^(25^^"^ - 2a2fo(ro - 3) - 2at;~^ - 3a^ + 2fo^ - 3fo^)5f^ (^rj, 



and 

/3(5) 



M'^5r 



|ro^ ro^5f ^ + Aofo^-^e^ + Ao (^^(fo + 2) + fo^) -^fA^ (8dfou'l'u\fo - 1) + lOa^u'u' 



SAo^fo^ 

+ a2[2fo(n'^)'(2-3fo)-5(n*)2 



(A9) 



5a 



{n^f - fo[roHu^? + 2{u')\fo - 2)])6f^ 



7.4:/~.(t>\2 



+ Ao(15a^(u' 



305Vu'^ + a2[fo(19fo - 6)(n'^)2 + 15(n*)2] + 2afouV(6 - llfo) 



+ fo[4fo^(n'^)2 + 3(u*)2(fo - 2)])(502 ^ Ao\*(3n* - 20^^50' 



2Ao^u*n^ 



52,i*(fo + 2) 



tn J- ^ 3 ~<^ 



2'u''a + fo n 



n* - att'^ 



3fo^<5r" + 3Aofo^<5e' + Aq a"[4foln'^)'(ro + 2) + 3fo' + 6fo 



Ui 



t\2] 



8u'd^{ro + 2)d'^ - Su'dfo^"^ + 25^(^0 + 2f{u'^)'^ + ro'^[2ro'^{u'^f + 3] 



*1- 

(AlO) 



We note that our expression for at^\ contains terms which are of 0{6x^) and therefore it is not 
strictly consistent with our original definition [a Taylor series truncated at 0{6x^)]. However, the 
inclusion of these terms allows us to express 0(5) more compactly, and we opt to retain them 
despite the slight abuse of notation. Of course, the additional 0{6x^) terms do not affect the 4th- 
order nature of the puncture. The 4th-order puncture employed in our numerical implementation 



uses 0(5) as it is given in Eq. (A9). A Mathematica workbook with details of the fourth-order 



puncture function is available [46jT 



Appendix B: Second-order effective source 

[21 [21 

For the 2nd-order effective source, defined by S^^ = S{x) — D^^ , we obtain the 7?i-mode 



decomposition 



sTit,r,e) = -e 



~im{nt+A(j)) ST^ a j 
k=l 



m 
k 1 



(Bl) 
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where A(/) was defined in Eq. (44) and the integrals /{", . . . , I5" are defined and expressed in terms of 
elliptic integrals in Eqs. (B9)-(B13) of Paper I. Since there is the possibility of notional confusion, 
let us give one of these integrals explicitly: 



rm 

-'1 = 



s-fe— • 



7 



M^B^^ 



[pTK{p)Ki7) + p-'pTE{p)E{7)], 



(B2) 



where S(3) is defined in Eq. (38), B is given in Eq. (45), and p and 7 are given in Eq. (47). The 
polynomials p^, p^ for m = 0, . . . ,5 are given in Table II of Ref. |25) : higher-order terms can be 
computed easily with a symbolic algebra package. 

The quantities Sk in (iBll) are Kerr-generalizations of the ones given in Paper I. They read 



Si = M'p''[{r-l)X + A{s2oo + 3s3ooSr) + {l + SOcote){so20 + Si2odr)], 



S2 
S3 

^5 



M'^p-^ 



T 



4^-2 



-3MV 



2nf^ + f7^7«j B-2{f- l)sio2 
y - S102 j + \Sq2Q + Si20'5r j 5(9^ 



4„-2 



3M*p 



As- 



102 



7" 



2fi7*'^ + ^2yt ) B 



-3M^p-^ 



X 



A ( y + S102 ) + ( S020 + si2o5f ] se"^ 



(B3) 
(B4) 

(B5) 

(B6) 

(B7) 



where the coefficients Sijk are the ones given in Appendix [Al a tilde denotes a-dimensionalization 
using M (e.g., f = r/M or (7 = Mfi), p^ is defined in (jol), and we have introduced 



X = 2s2oo5r + 3s3oo'^r^ + Si2o<50^ + 2s 



102, 



(B8) 



(r^ + g^)^ -2 . 2 
= a sm I 



7 
7' 
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2af 

"X' 

1 g^ 

sin^ 6* A ' 



(B9) 
(BIO) 
(Bll) 



Finally, S^ featuring in Eq. (24) is given by — rAp^S ^S'j.g . 
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